﻿module PCA
open MathNet.Numerics
open MathNet.Numerics.FSharp
open MathNet.Numerics.LinearAlgebra.Double
open System.Collections.Generic
open Utils
open Types


//Function that calculate the distance matrix between two matrix
//if X has M rows and N columns and C has L rows and N columns then
//the result matrix has dimension M x L
let dist (x : DenseMatrix) (c : DenseMatrix) =
  
  let a = (ones(c.RowCount) * sum_column ((matrix_square x).Transpose() :?> DenseMatrix)).Transpose()
  let b = ones(x.RowCount) * sum_column((matrix_square c).Transpose() :?> DenseMatrix)
  let z = 2.0 * x * (c.Transpose()) :?> DenseMatrix
  a + b - z :?> DenseMatrix

//Function that calculate the mean of a set 
let private mean (points : List<DenseVector>) =
  (points |> Seq.fold(fun sum point -> sum + point) (DenseVector(256))) / float points.Count

//function to compute the covariance of a list of points and a given mean
let compute_covariance_matrix (points : List<DenseVector>) (mean : DenseVector) =
  let res = 
    (points |> Seq.fold(fun sum point -> 
                          let diff = point - mean
                          let diff_trans = DenseVector.OuterProduct(diff, diff)
                          sum + diff_trans) (DenseMatrix(256)))
  res.Divide(float points.Count) :?> DenseMatrix

//compute the eigen vector
let eigen_vectors(matrix : DenseMatrix) =
  matrix.Svd(true).VT() :?> DenseMatrix

//function that compute the principal components vector
let alpha_matrix (pc_matrix : DenseMatrix) (mean : DenseVector) (training_points : List<DenseVector>) =
  let alpha_matrix = DenseMatrix(training_points.Count, pc_matrix.ColumnCount)
  for i in 0 .. alpha_matrix.RowCount - 1 do
    for j in 0 .. alpha_matrix.ColumnCount - 1 do
      alpha_matrix.[i,j] <-  (training_points.[i] - mean) * pc_matrix.Row(j)
  alpha_matrix

//compute the principal component alghoritm
let PCA (train : List<DenseVectorWrapper>) (test : List<DenseVectorWrapper>) = 
  let train = [for elem in train do yield elem.DenseVector]
  let test = [for elem in test do yield elem.DenseVector]
  let mean_train = mean (List(train))
  let s = compute_covariance_matrix (List(train)) mean_train
  let eigen_vectors = eigen_vectors s
  let alpha_train = alpha_matrix eigen_vectors mean_train (List(train))
  let alpha_test = alpha_matrix eigen_vectors mean_train (List(test))
  printf "PCA done!" 
  alpha_train, alpha_test

//compute the nearest neighbour alghoritm
let KNN (alpha_train  : DenseMatrix) 
        (alpha_test   : DenseMatrix) 
        (train_set    : List<DenseVectorWrapper>)
        (test_set     : List<DenseVectorWrapper>)
        (k : int) =

  //reduce the alpa train and test matrixes
  let beta_train = alpha_train.SubMatrix(0, alpha_train.RowCount, 0, k) :?> DenseMatrix
  let beta_test =  alpha_test.SubMatrix(0, alpha_test.RowCount, 0, k) :?> DenseMatrix

  //compute the distance matrix
  let dist_sq = dist beta_test beta_train 

  [for i in 0 .. test_set.Count - 1 do
    let r_i = dist_sq.Row(i)
    //takes the minimum of the distances
    let min_index = r_i.MinimumIndex()  
    //takes the value of the class corrisponding to the point of the training with 
    //the minimal distance  
    yield train_set.[min_index]]


//main test funtion
let main (get_train_and_test : unit -> List<DenseVectorWrapper> * List<DenseVectorWrapper>) =


  //number of cycles
  let cycles = 1

  let error_matrix = DenseMatrix(cycles,256)
  for j = 0 to cycles - 1 do

    do printfn "Cycle %d done!" j

    let train, test = get_train_and_test ()
    let alpha_train, alpha_test = PCA train test      
    //compute the knn for 255 times for 
    for i = 1 to 255 do
      let mutable total_error = 0
      let result = KNN alpha_train alpha_test train test i
      let result_list = Seq.zip test result
      for (data,classification) in result_list do
        if data.RealNumber <> classification.RealNumber then
          total_error <- total_error + 1
      error_matrix.[j,i] <- (float total_error) / (test |> Seq.length |> float)
    


  //compute the mean hit rate
  let error_vector =  DenseVector(256)
  for i = 1 to 255 do
    error_vector.[i] <- 1.0 - (error_matrix.Column(i) |> Seq.fold(fun s error -> s + error) 0.0) / (float (error_matrix.RowCount))

  //plot the results
  let input, alpha = 
    let fst,snd = get_train_and_test()
    let alpha_train, alpha_test = PCA fst snd
    fst.AddRange snd
    fst, concat_matrix alpha_train alpha_test
  do plot_with_data alpha input error_vector